Dielectric Breakdown of a Mott Insulator 
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We study the non-equilibrium steady state of a Mott insulator coupled to a thermostat and driven 
by a constant electric field, starting from weak fields, until the dielectric breakdown, and beyond. 
We find that the conventional Zener picture does not describe the steady-state physics. In particular, 
the current at weak field is found to be controlled by the dissipation. Moreover, in connection with 
the electric field driven dimensional crossover, we find that the dielectric breakdown occurs when 
the field strength is on the order of the Mott gap of the corresponding lower dimensional system. 
We also report a resonance and the melt-down of the quasi-particle peak when the field strength is 
half of this Mott gap. 



An early achievement in the understanding of the non- 
linear response of electronic systems driven by strong 
electric fields is due to Zener in 1934 [1|. He computed the 
rate of interband transitions of a one-dimensional non- 
interacting band insulator in a constant electric field, as- 
suming that there is no back-feeding from the conduction 
to the valence band. This predicted a threshold electric 
field Eth above which the dielectric breakdown of the in- 
sulator occurs. 

Following Oka, Arita, and Aoki's proposal that this 
single electron picture also applies to Mott insulators [2| , 
many efforts have been done to check their idea by test- 
ing Zener's formula: Eth °c A 2 where A is the gap of the 
insulator. Numerically, this out-of-equilibrium strongly 
interacting problem has been tackled by means of time- 
dependent (TD) methods such as TD density matrix 
renormalization group [3J or TD exact diagonalization 
in Id finite systems [4|, and by TD dynamical mean- 
field theory in infinite dimensions [5]. There, the lack of 
a dissipation mechanism (necessary to get a non-trivial 
steady-state as earlier understood by jy, Q), causes a 
continuous heating up of the system [8j . 

Experimentally, the electric field dependence of the 
current density is extracted from the current-voltage 
characteristic when applying a bias voltage on large 
samples [9j. Several examples exhibit a much smaller 
threshold field than the estimation from Zener's for- 
mula 10]. In this Letter, we address this problem by 
driving out of equilibrium a two-dimensional (2d) Hub- 
bard model coupled to a dissipative thermostat. We treat 
both the strong electric field and the strong interaction, 
and we bypass the transient dynamics by means of the 
non-equilibrium steady-state dynamical mean-field the- 
ory (NESS-DMFT) developed recently by the author and 
collaborators fllj . 

Hereafter, we describe the model and detail the com- 
putations. Then, we summarize the influence of the dis- 
sipation on the equilibrium physics of the Mott transi- 
tion. Later, we study the influence of the electric field 
on the spectral properties of the Mott insulator and ar- 
gue that the dissipation is the leading mechanism for the 
interband current. Afterwards, we undertake the system- 



atic exploration of the non-linear response of the system 
as the electric field is increased and as the dimensional 
crossover to the corresponding Id system takes place, un- 
til the full dimensional reduction predicted on general 
grounds in [11]. In particular, we discuss a small jump in 
the conductivity and the melt-down of the quasi-particle 
peak when the field strength is half of the Mott gap of 
this Id system. We also detail the physics of the dielec- 
tric breakdown that is found when the field strength is on 
the order of this Mott gap, contrary to Zener's picture. 

Model. We consider the Hubbard model on a d = 2 
square lattice. The static and uniform electric field is 
set along an axis of the lattice: E = Eu x with E > 0. 
The Lagrangian of the system coupled to its environment 
reads (we set h = 1 and use the conventions of [llj) 
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where a a and Ci a are the Grassmann fields representing 
an electron at site i with spin a G {j, !}. U is the on-site 
Coulombic interaction and iy = (a/27r) 2 Jdk e lk ' XiJ e(k) 
sets the hopping amplitude between two nearest neigh- 
bors distant of a: e(k) = eo [cos{k x a) + cos(k y a)], each 
dimension contributing by 2eo to the bandwidth of the 
equilibrium non-interacting (E = U = 0) system. Inte- 
grals over k x and k y are computed between —ir/a and 
tk I a. The last term in ([I) is the coupling to the thermo- 
stat which is composed of independent non-interacting 
electronic reservoirs in equilibrium at a very low temper- 
ature T and a chemical potential fj,o — U/2 in order to 
work at half-filling i.e. with one electron per site in av- 
erage (we also restrict ourselves to the paramagnetic so- 
lution and drop the spin indices) . 7 is a real hopping pa- 
rameter, the 6's represent the electrons in the reservoirs, 
and I labels their energy levels. The Peierls phase factors, 
(Xij(t) = qf^dx ■ A(i,x) and 0i(t) = / dt' (j) l {t') 1 are 
required by the gauged U(l) symmetry associated with 
the conservation of the charge q of the electrons. <j> and 



A are the scalar and vector potentials: E = — V</> — 9* A. 
\q\Ea is the energy an electron acquires when hopping 
to a neighboring site under the work of the electric field. 
To work with gauge-invariant quantities, we use the vari- 
ables vj = lu — (j) and k = k — qA. and later absorb the 
Hartree shift by redefining w — U/2 into w. 

An efficient dissipation is achieved if the bandwidth 
W of the local density of states (DOS) of the reservoirs 
is the largest energy scale. The other details of these 
DOS are not relevant and we take them to be Gaussian, 
yielding a contribution of the dissipation to the Keldysh 
self-energy: £f h (za) = T exp(-w 2 /nW 2 ) tanh(tn/2fc B T) 
where V = j 2 /W. We work at small dissipation T 
but large enough for the momentum resolved spectral 
function to be positive everywhere, ensuring a stable 
steady state. Otherwise, this signals oscillatory insta- 
bilities (such as Bloch oscillations) developing on top of 



12] 



the steady-state solution 

Computational details. The non-equilibrium steady 
state is solved by means of the NESS-DMFT algo- 
rithm developed in [ll| and based on a gauge-invariant 
Schwinger-Kcldysh formalism [l3| |. The interaction con- 
tribution to the retarded and Keldysh self-energies (£^ 
and Ey ) are computed using second order iterated per- 
turbation theory in U (IPT) as the impurity solver. Al- 
though it is not a ^-derivable approximation, it is a cur- 
rent conserving approximation at half-filling 14]. For 
each value of the electric field, the dressed retarded and 
Keldysh Green's functions (G§ and G^ ) are obtained in 
the strongly interacting regime by starting from the non- 
interacting solution, then by slowly increasing the inter- 
action ([/ H> U + 5U) while converging at each step the 
impurity and the following lattice equations [see Eqs. (7), 
(8), (9), and (3) in [Hi], 
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where 5T, 

take further advantage of both the mean-field approxi- 
mation and the geometry of the setup, the star product 
is evaluated in the mixed (tn; n x ; Kj / )-space where n x G Z: 

[/ * g] (zu; n x ; n y ) =^ / {w + m x qEa/2; n x - m x ; K y ) 

x g(w + (m x - n x )qEa/2;m x ;n y ) , (4) 

with f(m;n x ;n y ) = (a/2n) Jdn x & Kxn * a f (vo , k). Each 
evaluation of Eqs. (2) and (3) requires performing a sin- 
gle numerical summation and the overall computation 
is slower than the equilibrium algorithm by a factor 
2N X = 2 x 27r / ' a8n x only, where Sk x is the discretization 
step for k x . Hereafter, numerical results are obtained 
with eg = a = q = Icb = I- 

Influence of the dissipation in equilibrium. In equilib- 
rium {E — 0) and as the interaction U is increased, the 
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FIG. 1: In-gap DOS for U = 20 and E = 6 (T = 0.05). Dotted 
line: very weak dissipation T = 0.09 (case presenting oscillatory 
instabilities) revealing the BZ "archipelagos" centered around e a = 
±qEa and both composed of two Id BZ islands at e a ± U/2. Solid 
line: the same for T = 0.25 where the quasi-particle peak around 
e = is much stronger and the details of the BZ islands are now 
almost indistinguishable. 



Hubbard model exhibits a well-known quantum phase 
transition from a metal to a Mott insulator characterized 
by the opening of an energy gap A ~ U — 2den separating 
the so-called Hubbard bands (l5( . The presence of a weak 
dissipation T smoothens the sharp features of the spec- 
tral function over an energy window T. In particular, the 
edges of the Hubbard bands leak into the gap, respon- 
sible for a dissipative in-gap DOS controlled by T/U 2 . 
Dissipation also delays the transition which turns into a 
smooth crossover taking place onto an extended region 
in U [19( . There, what is left of the metal manifests itself 
by a weakly dispersive Kondo-likc resonance of width cok, 
centered around the Fermi level, and containing a frac- 
tion Z of all the states. When increasing U, the height 
of this quasi-particle peak is first roughly constant (and 
decreases with T) while u>k decreases continuously. Deep 
in the strongly interacting phase, the peak becomes con- 
trolled by the dissipation as lok is rather constant (and 
set by r) while its height vanishes as l/U 2 (and grows 
with r) [see its dependence on U and T in Fig. [21(a)]. 

Influence of the electric field on the spectral proper- 
ties. Since deep in the strongly interacting regime, each 
Hubbard band exhibits the spectral features of a single 
non-interacting band and U only enters through the gap 
A [see Fig.[3Ja)], one can expect the effects of the electric 
field on the spectral function of the Mott insulator to be 
similar to the case of a non-interacting band insulator. 
In the well-known case of a single non-interacting band, 
some Bloch-Zener (BZ) islands appear in the DOS be- 
yond the edges of the band, equally spaced in energy by 
l^l-Ba, with a weight that is exponentially killed on a scale 
e (\q\Ea) 2 / 3 as one gets away from the band edges, and 
the energy structure of which is controlled by the DOS 
of the equilibrium Id system along the y-direction 12j. 
In our Mott insulator case, we find that this scenario in- 
deed occurs as we observe similar islands in the DOS. 
They have the structure of the corresponding equilib- 
rium Id Mott insulator. Since the latter has a gapped 
DOS, the islands are in fact "archipelagos" centered on 
multiples of ±qEa and composed of two islands of width 
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FIG. 2: (color online) (a) Height of the equilibrium quasi-particlc 
peak, p(e = 0), as a function of U for different dissipations F (E = 
0, T = 0.05). (b) The two melt-downs of the quasi-particle peak 
at E < ojjc and E ~ U/2, followed by the growth of the peak 
of the equilibrium Id model. Dashed line: level of the dissipativc 
background at e ~ estimated from the equilibrium data (U = 
16, T= 0.05, r = 0.25). 



2eo and separated by U. We illustrate in Fig. Q] this 
rich structure of the DOS between the Hubbard bands. 
These in-gap islands allow the transition of carriers from 
the lower to the upper Hubbard band by successive exci- 
tations of energy |g|.Ea. However, the dissipation creates 
a continuous in-gap DOS, damped as a power law as one 
gets away from the band edges, and it is therefore ex- 
pected to be the main contribution for those in-gap states 
(see Fig.[T]). It was indeed the case in all the stable steady 
states we expolred. 

Before we start the systematic study of the non-linear 
regime, notice that as the field is increased, the system 
experiences a dimensional crossover from the insulating 
phase of the 2d equilibrium Hubbard model (at E = 0) to 
the insulating phase of the Id equilibrium model (when 
\q\Ea is the largest energy scale) [ll|. Since both exhibit 
a similar DOS deep in the strongly interacting regime (at 
least for the paramagnetic solutions obtained with the lo- 
cal approximation of the single-site DMFT), one expects 
a smooth variation from the 2d DOS with a pseudo gap 
A2d — U — Aeo between bands of width 4eo towards the Id 
DOS with a pseudo gap Aid — U — 2eo between bands of 
width 2eo |17| . Therefore, the qualitative features of the 
current characteristic can already be predicted by simply 
reasoning on Fig. [TJ 

Below, we detail the fate of the insulating phase when 
increasing the electric field by focusing on the momentum 
resolved spectral function p(e, k), the local DOS p(e), and 
the current density J(E) plotted in Fig. [5] The latter 
also provides qualitative informations on the asymmetry 
in k x of the momentum distribution function n(n): J ex 
— 2 fdn d K , x e(K) n(/c). 

\q\Ea <C 2e <C U. Let us start with very weak fields. 
At moderate values of U for which the quasi-particle peak 
is still present (Z > 0), a small electric field such as 
l^l-Ea < cok can excite the states lying in an energy shell 
ujk below the Fermi level (e = 0) to the empty states 
above, resulting in a tiny current. This reorganization 
of the distribution of occupied states around e = is 
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FIG. 3: (color online) (a) Equilibrium spectral function p(e, re) 
integrated over K y for U = 20. (b) The same for 2eo = 2 < E = 
8 < U = 20 is now almost re^-independent and the bands have a 
width 2e , similarly to the Id model. (T = 0.05, T = 0.35). 



qualitatively similar to having an effective temperature, 
and causes the partial melt-down of the quasi-particle 
peak [see Fig. [U(b)] 16]. As soon as \q\Ea is larger than 
lok, the transition rate is now controlled by the small in- 
gap DOS created by the leakage of the Hubbard bands, 
leading to a drop in the differential conductivity. Deep in 
the strongly interacting regime, the quasi-particle peak 
vanishes (Z —> 0) and the growth of the tiny current is 
controlled by the dissipative in-gap DOS which is on the 
order of T/A 2 . 

2eo <C \q\Ea « U. As the electric field intensity 
gets larger 2eo (i.e. the fraction of the non-interacting 
bandwidth corresponding to the x direction), p(e, k) loses 
much of its dependence on k x and becomes essentially the 
one of the Id Hubbard model [see Fig.^b)] [Sj. This is 
a first step towards the full dimensional reduction of the 
system. Meanwhile, n(n) is still very close to the one of 
the 2d Hubbard model in equilibrium [see Fig. HJa)] and 
the current is very weak. 

2en <C \q\Ea ~ U/2. When the electric field intensity 
is comparable with the energy separating the lower Hub- 
bard band with the Fermi level (e = 0), \q\Ea ~ Aid/2 ~ 
U/2 — eo, carriers can be excited from the former to 
the dissipative background around the latter. Concomi- 
tantly, the large amount of vacant states offered by the 
upper band favors a rapid pumping of these newly occu- 
pied states to the upper band. These resonant processes 
contribute to a significant increase of the current density 
until \q\Ea ~ U/2. Here again, a stronger dissipation fa- 
vors a larger current via the increase of the in-gap DOS. 
Notice also that the combination of the BZ archipela- 
gos centered at ±qEa creates a large mid-gap BZ island 
on top of the dissipative background (see Fig. [T]) which 
also contributes to this resonance. Furthermore, the re- 
organization of the distribution of occupied states around 
e = (the fraction of occupied states decreases signifi- 
cantly just below e = while it increases symmetrically 
above) is qualitatively similar to having a high effective 
temperature. This explains the complete melt-down of 
the quasi-particle peak that we observe until the reso- 
nance is broken when \q\Ea > U/2 + eo [see Fig.[2jb)]. 

2eo <C \q\Ea ~ U. When the electric field is of the 
magnitude of the Mott gap, \q\Ea ~ Aid, carriers in the 
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FIG. 4: (color online) (a) The momentum distribution function 
n(n) for 2eo = 2<_E = 8<C7 = 20is similar to 2d equilibrium 
[contrary to p(e, k,) in Fig. IH^b)] . (b) n(fc) just after the dielectric 
breakdown for 2e = 2 < E = 20 ~ (7 = 20. (T = 0.05, T = 0.35). 
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FIG. 5: (color online) (a) Current density J(E) for different U 
(T = 0.05, T = 0.25). The first jump is located at E ~ (7/2 and 
the maximum at E ~ £7. (b) The same data is plotted against 
E — U to prove the scaling J(E — (7) in the metalized regime. The 
dotted curve corresponds to T = 0.50. 



lower band can directly populate the upper band. This is 
the dielectric breakdown. The electric current increases 
rapidly as the field is further increased and reaches a 
maximum at \q\Ea ~ U. After the dielectric break- 
down, one expects U to be quite irrelevant and the scaling 
,/ = J(E — U) to hold [see Fig. Etb)] since the structure 
and the occupation of each Hubbard band are almost 
independent of U, except for small corrections due to 
the dissipative background and the quasi-particle peak 
(if any). The Hubbard bands can be seen as two non- 
interacting systems connected to a thermostat: electrons 
are excited from the first one to the second, then are ab- 
sorbed by the thermostat which also repopulates the first 
system. The dissipation enters this picture in two ways. 
One the one hand, a stronger dissipation accelerates the 
repopulatation of the lower Hubbard band and should 
therefore favor a larger current. On the other hand, the 
dissipation is expected to reduce the current because it is 
responsible for fewer states in the Hubbard bands since 
they leak into the gap and since it also strengthens the 
quasi-particle peak. All together, we show in Fig. EJb) 
that a stronger dissipation favors a smaller value of the 
maximum current but a larger current away from this 
maximum. Together with the sharp current increase, 
the weight of n(ft) is strongly displaced along k x [see 
Fig. HJb)]. Notice the sharper discontinuity of n(«) due 
to the fact that the lower (upper) Hubbard band has now 
a sizable fraction of unoccupied (occupied) states. 

2eo C(7< \q\Ea. When the electric field is stronger 
than any other energy scale, the dimensional reduction 
predicts that the system behaves as a collection of uncou- 
pled Id Hubbard chains in equilibrium ll(. The DOS 
being bounded, the electric field is too strong for any 
transition to take place (except in the outer dissipative 
background) as soon as \q\Ea > U + 2e - Both p(e, k) 
and n(K) are the ones of the Id Hubbard model in equi- 
librium and the current vanishes accordingly. 

Discussion. We have investigated the steady-state 
physics of a 2d Mott insulator driven out of equilibrium 
by a constant electric field and coupled to a thermostat. 
We argued that the interband current is mostly due to the 
presence of in-gap states created by the dissipation. Also 



contrary to Zener's picture, we observed the dielectric 
breakdown of the Mott insulator after \q\Ea ~ Ai^. Fur- 
thermore, we revealed a resonance around |q r |.E'a ~ U/2 
responsible for a small jump in the conductivity and for 
the melt-down of the quasi-particle peak. We also showed 
that the dimensional crossover takes place on two sepa- 
rated energy scales: the spectral properties turn to the 
ones of the Id Mott insulator as soon as \q\Ea 3> 2eo, 
whereas the distribution functions only reach thermal 
equilibrium in Id when |g|Sa 3> U . 

We expect this scenario to be also relevant for 3d sam- 
ples crossing over to 2d, where the DMFT solutions are 
all the more valid. We also believe that our results can 
be put to experimental test with cold atoms trapped in 
optical lattices where strong electric fields (|q r |.E'a > U) 
can be mimicked by forcing the lattice potential [18| and 
the dissipation can be engineered by coupling the Mott 
insulator to a superfluid fraction of the atomic conden- 
sate. 
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